home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 16 / CU Amiga Magazine's Super CD-ROM 16 (1997-10-16)(EMAP Images)(GB)[!][issue 1997-11].iso / CUCD / Graphics / Ghostscript / source / gsmatrix.c < prev    next >
C/C++ Source or Header  |  1996-12-07  |  13KB  |  427 lines
  1. /* Copyright (C) 1989, 1995, 1996 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of Aladdin Ghostscript.
  4.   
  5.   Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  6.   or distributor accepts any responsibility for the consequences of using it,
  7.   or for whether it serves any particular purpose or works at all, unless he
  8.   or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  9.   License (the "License") for full details.
  10.   
  11.   Every copy of Aladdin Ghostscript must include a copy of the License,
  12.   normally in a plain ASCII text file named PUBLIC.  The License grants you
  13.   the right to copy, modify and redistribute Aladdin Ghostscript, but only
  14.   under certain conditions described in the License.  Among other things, the
  15.   License requires that the copyright notice and this notice be preserved on
  16.   all copies.
  17. */
  18.  
  19. /* gsmatrix.c */
  20. /* Matrix operators for Ghostscript library */
  21. #include "math_.h"
  22. #include "gx.h"
  23. #include "gserrors.h"
  24. #include "gxfarith.h"
  25. #include "gxfixed.h"
  26. #include "gxmatrix.h"
  27.  
  28. /* The identity matrix */
  29. private gs_matrix gs_identity_matrix =
  30.     { identity_matrix_body };
  31.  
  32. /* ------ Matrix creation ------ */
  33.  
  34. /* Create an identity matrix */
  35. void
  36. gs_make_identity(gs_matrix *pmat)
  37. {    *pmat = gs_identity_matrix;
  38. }
  39.  
  40. /* Create a translation matrix */
  41. int
  42. gs_make_translation(floatp dx, floatp dy, gs_matrix *pmat)
  43. {    *pmat = gs_identity_matrix;
  44.     pmat->tx = dx;
  45.     pmat->ty = dy;
  46.     return 0;
  47. }
  48.  
  49. /* Create a scaling matrix */
  50. int
  51. gs_make_scaling(floatp sx, floatp sy, gs_matrix *pmat)
  52. {    *pmat = gs_identity_matrix;
  53.     pmat->xx = sx;
  54.     pmat->yy = sy;
  55.     return 0;
  56. }
  57.  
  58. /* Create a rotation matrix. */
  59. /* The angle is in degrees. */
  60. int
  61. gs_make_rotation(floatp ang, gs_matrix *pmat)
  62. {    gs_sincos_t sincos;
  63.     gs_sincos_degrees(ang, &sincos);
  64.     pmat->yy = pmat->xx = sincos.cos;
  65.     pmat->xy = sincos.sin;
  66.     pmat->yx = -sincos.sin;
  67.     pmat->tx = pmat->ty = 0.0;
  68.     return 0;
  69. }
  70.  
  71. /* ------ Matrix arithmetic ------ */
  72.  
  73. /* Multiply two matrices.  We should check for floating exceptions, */
  74. /* but for the moment it's just too awkward. */
  75. /* Since this is used heavily, we check for shortcuts. */
  76. int
  77. gs_matrix_multiply(const gs_matrix *pm1, const gs_matrix *pm2, gs_matrix *pmr)
  78. {    double xx1 = pm1->xx, yy1 = pm1->yy;
  79.     double tx1 = pm1->tx, ty1 = pm1->ty;
  80.     double xx2 = pm2->xx, yy2 = pm2->yy;
  81.     double xy2 = pm2->xy, yx2 = pm2->yx;
  82.     if ( is_xxyy(pm1) )
  83.     {    pmr->tx = tx1 * xx2 + pm2->tx;
  84.         pmr->ty = ty1 * yy2 + pm2->ty;
  85.         if ( is_fzero(xy2) )
  86.           pmr->xy = 0;
  87.         else
  88.           pmr->xy = xx1 * xy2,
  89.           pmr->ty += tx1 * xy2;
  90.         pmr->xx = xx1 * xx2;
  91.         if ( is_fzero(yx2) )
  92.           pmr->yx = 0;
  93.         else
  94.           pmr->yx = yy1 * yx2,
  95.           pmr->tx += ty1 * yx2;
  96.         pmr->yy = yy1 * yy2;
  97.     }
  98.     else
  99.     {    double xy1 = pm1->xy, yx1 = pm1->yx;
  100.         pmr->xx = xx1 * xx2 + xy1 * yx2;
  101.         pmr->xy = xx1 * xy2 + xy1 * yy2;
  102.         pmr->yy = yx1 * xy2 + yy1 * yy2;
  103.         pmr->yx = yx1 * xx2 + yy1 * yx2;
  104.         pmr->tx = tx1 * xx2 + ty1 * yx2 + pm2->tx;
  105.         pmr->ty = tx1 * xy2 + ty1 * yy2 + pm2->ty;
  106.     }
  107.     return 0;
  108. }
  109.  
  110. /* Invert a matrix.  Return gs_error_undefinedresult if not invertible. */
  111. int
  112. gs_matrix_invert(const gs_matrix *pm, gs_matrix *pmr)
  113. {    /* We have to be careful about fetch/store order, */
  114.     /* because pm might be the same as pmr. */
  115.     if ( is_xxyy(pm) )
  116.     {    if ( is_fzero(pm->xx) || is_fzero(pm->yy) )
  117.           return_error(gs_error_undefinedresult);
  118.         pmr->tx = - (pmr->xx = 1.0 / pm->xx) * pm->tx;
  119.         pmr->xy = 0.0;
  120.         pmr->yx = 0.0;
  121.         pmr->ty = - (pmr->yy = 1.0 / pm->yy) * pm->ty;
  122.     }
  123.     else
  124.     {    double det = pm->xx * pm->yy - pm->xy * pm->yx;
  125.         double mxx = pm->xx, mtx = pm->tx;
  126.         if ( det == 0 )
  127.           return_error(gs_error_undefinedresult);
  128.         pmr->xx = pm->yy / det;
  129.         pmr->xy = - pm->xy / det;
  130.         pmr->yx = - pm->yx / det;
  131.         pmr->yy = mxx / det;    /* xx is already changed */
  132.         pmr->tx = - (mtx * pmr->xx + pm->ty * pmr->yx);
  133.         pmr->ty = - (mtx * pmr->xy + pm->ty * pmr->yy);    /* tx ditto */
  134.     }
  135.     return 0;
  136. }
  137.  
  138. /* Translate a matrix, possibly in place. */
  139. int
  140. gs_matrix_translate(const gs_matrix *pm, floatp dx, floatp dy, gs_matrix *pmr)
  141. {    gs_point trans;
  142.     int code = gs_distance_transform(dx, dy, pm, &trans);
  143.     if ( code < 0 )
  144.       return code;
  145.     if ( pmr != pm )
  146.       *pmr = *pm;
  147.     pmr->tx += trans.x;
  148.     pmr->ty += trans.y;
  149.     return 0;
  150. }
  151.  
  152. /* Scale a matrix, possibly in place. */
  153. int
  154. gs_matrix_scale(const gs_matrix *pm, floatp sx, floatp sy, gs_matrix *pmr)
  155. {    pmr->xx = pm->xx * sx;
  156.     pmr->xy = pm->xy * sx;
  157.     pmr->yx = pm->yx * sy;
  158.     pmr->yy = pm->yy * sy;
  159.     if ( pmr != pm )
  160.       { pmr->tx = pm->tx;
  161.         pmr->ty = pm->ty;
  162.       }
  163.     return 0;
  164. }
  165.  
  166. /* Rotate a matrix, possibly in place.  The angle is in degrees. */
  167. int
  168. gs_matrix_rotate(const gs_matrix *pm, floatp ang, gs_matrix *pmr)
  169. {    double mxx, mxy;
  170.     gs_sincos_t sincos;
  171.     gs_sincos_degrees(ang, &sincos);
  172.     mxx = pm->xx, mxy = pm->xy;
  173.     pmr->xx = sincos.cos * mxx + sincos.sin * pm->yx;
  174.     pmr->xy = sincos.cos * mxy + sincos.sin * pm->yy;
  175.     pmr->yx = sincos.cos * pm->yx - sincos.sin * mxx;
  176.     pmr->yy = sincos.cos * pm->yy - sincos.sin * mxy;
  177.     if ( pmr != pm )
  178.       { pmr->tx = pm->tx;
  179.         pmr->ty = pm->ty;
  180.       }
  181.     return 0;
  182. }
  183.  
  184. /* ------ Coordinate transformations (floating point) ------ */
  185.  
  186. /* Note that all the transformation routines take separate */
  187. /* x and y arguments, but return their result in a point. */
  188.  
  189. /* Transform a point. */
  190. int
  191. gs_point_transform(floatp x, floatp y, const gs_matrix *pmat,
  192.   gs_point *ppt)
  193. {    ppt->x = x * pmat->xx + pmat->tx;
  194.     ppt->y = y * pmat->yy + pmat->ty;
  195.     if ( !is_fzero(pmat->yx) )
  196.       ppt->x += y * pmat->yx;
  197.     if ( !is_fzero(pmat->xy) )
  198.       ppt->y += x * pmat->xy;
  199.     return 0;
  200. }
  201.  
  202. /* Inverse-transform a point. */
  203. /* Return gs_error_undefinedresult if the matrix is not invertible. */
  204. int
  205. gs_point_transform_inverse(floatp x, floatp y, const gs_matrix *pmat,
  206.   gs_point *ppt)
  207. {    if ( is_xxyy(pmat) )
  208.        {    if ( is_fzero(pmat->xx) || is_fzero(pmat->yy) )
  209.           return_error(gs_error_undefinedresult);
  210.         ppt->x = (x - pmat->tx) / pmat->xx;
  211.         ppt->y = (y - pmat->ty) / pmat->yy;
  212.         return 0;
  213.        }
  214.     else if ( is_xyyx(pmat) )
  215.        {    if ( is_fzero(pmat->xy) || is_fzero(pmat->yx) )
  216.           return_error(gs_error_undefinedresult);
  217.         ppt->x = (y - pmat->ty) / pmat->xy;
  218.         ppt->y = (x - pmat->tx) / pmat->yx;
  219.         return 0;
  220.        }
  221.     else
  222.        {    /* There are faster ways to do this, */
  223.         /* but we won't implement one unless we have to. */
  224.         gs_matrix imat;
  225.         int code = gs_matrix_invert(pmat, &imat);
  226.         if ( code < 0 )
  227.           return code;
  228.         return gs_point_transform(x, y, &imat, ppt);
  229.        }
  230. }
  231.  
  232. /* Transform a distance. */
  233. int
  234. gs_distance_transform(floatp dx, floatp dy, const gs_matrix *pmat,
  235.   gs_point *pdpt)
  236. {    pdpt->x = dx * pmat->xx;
  237.     pdpt->y = dy * pmat->yy;
  238.     if ( !is_fzero(pmat->yx) )
  239.       pdpt->x += dy * pmat->yx;
  240.     if ( !is_fzero(pmat->xy) )
  241.       pdpt->y += dx * pmat->xy;
  242.     return 0;
  243. }
  244.  
  245. /* Inverse-transform a distance. */
  246. /* Return gs_error_undefinedresult if the matrix is not invertible. */
  247. int
  248. gs_distance_transform_inverse(floatp dx, floatp dy,
  249.   const gs_matrix *pmat, gs_point *pdpt)
  250. {    if ( is_xxyy(pmat) )
  251.        {    if ( is_fzero(pmat->xx) || is_fzero(pmat->yy) )
  252.           return_error(gs_error_undefinedresult);
  253.         pdpt->x = dx / pmat->xx;
  254.         pdpt->y = dy / pmat->yy;
  255.        }
  256.     else if ( is_xyyx(pmat) )
  257.       {    if ( is_fzero(pmat->xy) || is_fzero(pmat->yx) )
  258.           return_error(gs_error_undefinedresult);
  259.         pdpt->x = dy / pmat->xy;
  260.         pdpt->y = dx / pmat->yx;
  261.       }
  262.     else
  263.        {    double det = pmat->xx * pmat->yy - pmat->xy * pmat->yx;
  264.         if ( det == 0 )
  265.           return_error(gs_error_undefinedresult);
  266.         pdpt->x = (dx * pmat->yy - dy * pmat->yx) / det;
  267.         pdpt->y = (dy * pmat->xx - dx * pmat->xy) / det;
  268.        }
  269.     return 0;
  270. }
  271.  
  272. /* Compute the bounding box of 4 points. */
  273. int
  274. gs_points_bbox(const gs_point pts[4], gs_rect *pbox)
  275. {
  276. #define assign_min_max(vmin, vmax, v0, v1)\
  277.   if ( v0 < v1 ) vmin = v0, vmax = v1; else vmin = v1, vmax = v0
  278. #define assign_min_max_4(vmin, vmax, v0, v1, v2, v3)\
  279.   { double min01, max01, min23, max23;\
  280.     assign_min_max(min01, max01, v0, v1);\
  281.     assign_min_max(min23, max23, v2, v3);\
  282.     vmin = min(min01, min23);\
  283.     vmax = max(max01, max23);\
  284.   }
  285.     assign_min_max_4(pbox->p.x, pbox->q.x,
  286.              pts[0].x, pts[1].x, pts[2].x, pts[3].x);
  287.     assign_min_max_4(pbox->p.y, pbox->q.y,
  288.              pts[0].y, pts[1].y, pts[2].y, pts[3].y);
  289. #undef assign_min_max
  290. #undef assign_min_max_4
  291.     return 0;
  292. }
  293.  
  294. /* Transform or inverse-transform a bounding box. */
  295. /* Return gs_error_undefinedresult if the matrix is not invertible. */
  296. private int
  297. bbox_transform_either_only(const gs_rect *pbox_in, const gs_matrix *pmat,
  298.   gs_point pts[4],
  299.   int (*point_xform)(P4(floatp, floatp, const gs_matrix *, gs_point *)))
  300. {    int code;
  301.     if ( (code = (*point_xform)(pbox_in->p.x, pbox_in->p.y, pmat, &pts[0])) < 0 ||
  302.          (code = (*point_xform)(pbox_in->p.x, pbox_in->q.y, pmat, &pts[1])) < 0 ||
  303.          (code = (*point_xform)(pbox_in->q.x, pbox_in->p.y, pmat, &pts[2])) < 0 ||
  304.          (code = (*point_xform)(pbox_in->q.x, pbox_in->q.y, pmat, &pts[3])) < 0
  305.        )
  306.       DO_NOTHING;
  307.     return code;
  308. }
  309.  
  310. private int
  311. bbox_transform_either(const gs_rect *pbox_in, const gs_matrix *pmat,
  312.   gs_rect *pbox_out,
  313.   int (*point_xform)(P4(floatp, floatp, const gs_matrix *, gs_point *)))
  314. {    int code;
  315.     /*
  316.      * In principle, we could transform only one point and two
  317.      * distance vectors; however, because of rounding, we will only
  318.      * get fully consistent results if we transform all 4 points.
  319.      * We must compute the max and min after transforming,
  320.      * since a rotation may be involved.
  321.      */
  322.     gs_point pts[4];
  323.     if ( (code = bbox_transform_either_only(pbox_in, pmat, pts, point_xform)) < 0 )
  324.       return code;
  325.     return gs_points_bbox(pts, pbox_out);
  326. }
  327. int
  328. gs_bbox_transform(const gs_rect *pbox_in, const gs_matrix *pmat,
  329.   gs_rect *pbox_out)
  330. {    return bbox_transform_either(pbox_in, pmat, pbox_out,
  331.                      gs_point_transform);
  332. }
  333. int
  334. gs_bbox_transform_only(const gs_rect *pbox_in, const gs_matrix *pmat,
  335.   gs_point points[4])
  336. {    return bbox_transform_either_only(pbox_in, pmat, points,
  337.                       gs_point_transform);
  338. }
  339. int
  340. gs_bbox_transform_inverse(const gs_rect *pbox_in, const gs_matrix *pmat,
  341.   gs_rect *pbox_out)
  342. {    return bbox_transform_either(pbox_in, pmat, pbox_out,
  343.                      gs_point_transform_inverse);
  344. }
  345.  
  346. /* ------ Coordinate transformations (to fixed point) ------ */
  347.  
  348. #define f_fits_in_fixed(f) f_fits_in_bits(f, fixed_int_bits)
  349.  
  350. /* Transform a point with a fixed-point result. */
  351. int
  352. gs_point_transform2fixed(const gs_matrix_fixed *pmat,
  353.   floatp x, floatp y, gs_fixed_point *ppt)
  354. {    fixed px, py, t;
  355.     double dtemp;
  356.     int code;
  357.     if ( !pmat->txy_fixed_valid )
  358.       {    /* The translation is out of range.  Do the */
  359.         /* computation in floating point, and convert to */
  360.         /* fixed at the end. */
  361.         gs_point fpt;
  362.         gs_point_transform(x, y, (const gs_matrix *)pmat, &fpt);
  363.         if ( !(f_fits_in_fixed(fpt.x) && f_fits_in_fixed(fpt.y)) )
  364.           return_error(gs_error_limitcheck);
  365.         ppt->x = float2fixed(fpt.x);
  366.         ppt->y = float2fixed(fpt.y);
  367.         return 0;
  368.       }
  369.     if ( !is_fzero(pmat->xy) )
  370.       {    /* Hope for 90 degree rotation */
  371.         if ( (code = set_dfmul2fixed_vars(px, y, pmat->yx, dtemp)) < 0 ||
  372.              (code = set_dfmul2fixed_vars(py, x, pmat->xy, dtemp)) < 0
  373.            )
  374.           return code;
  375.         if ( !is_fzero(pmat->xx) )
  376.           { if ( (code = set_dfmul2fixed_vars(t, x, pmat->xx, dtemp)) < 0 )
  377.               return code;
  378.             px += t;        /* should check for overflow */
  379.           }
  380.         if ( !is_fzero(pmat->yy) )
  381.           { if ( (code = set_dfmul2fixed_vars(t, y, pmat->yy, dtemp)) < 0 )
  382.               return code;
  383.             py += t;        /* should check for overflow */
  384.           }
  385.       }
  386.     else
  387.       {    if ( (code = set_dfmul2fixed_vars(px, x, pmat->xx, dtemp)) < 0 ||
  388.              (code = set_dfmul2fixed_vars(py, y, pmat->yy, dtemp)) < 0
  389.            )
  390.           return code;
  391.         if ( !is_fzero(pmat->yx) )
  392.           { if ( (code = set_dfmul2fixed_vars(t, y, pmat->yx, dtemp)) < 0 )
  393.               return code;
  394.             px += t;        /* should check for overflow */
  395.           }
  396.       }
  397.     ppt->x = px + pmat->tx_fixed;    /* should check for overflow */
  398.     ppt->y = py + pmat->ty_fixed;    /* should check for overflow */
  399.     return 0;
  400. }
  401.  
  402. /* Transform a distance with a fixed-point result. */
  403. int
  404. gs_distance_transform2fixed(const gs_matrix_fixed *pmat,
  405.   floatp dx, floatp dy, gs_fixed_point *ppt)
  406. {    fixed px, py, t;
  407.     double dtemp;
  408.     int code;
  409.     if ( (code = set_dfmul2fixed_vars(px, dx, pmat->xx, dtemp)) < 0 ||
  410.          (code = set_dfmul2fixed_vars(py, dy, pmat->yy, dtemp)) < 0
  411.        )
  412.       return code;
  413.     if ( !is_fzero(pmat->yx) )
  414.       {    if ( (code = set_dfmul2fixed_vars(t, dy, pmat->yx, dtemp)) < 0 )
  415.           return code;
  416.         px += t;        /* should check for overflow */
  417.       }
  418.     if ( !is_fzero(pmat->xy) )
  419.       {    if ( (code = set_dfmul2fixed_vars(t, dx, pmat->xy, dtemp)) < 0 )
  420.           return code;
  421.         py += t;        /* should check for overflow */
  422.       }
  423.     ppt->x = px;
  424.     ppt->y = py;
  425.     return 0;
  426. }
  427.